
# Define an AEP class
setClass("AEP", representation(
	data = "data.frame",
	pobs_data = "matrix",
	nb_obs = "numeric",
	nb_var = "numeric",
	EmpiricalCDF = "list",
	ECDF = "function",
	N = "numeric",	
	binary_matrix = "array",
	sum_ones = "array",
	Hmultiplier = "array",
	alpha = "numeric",
	nb_iteration = "numeric",
	m = "array"
	)
)
setMethod("initialize",
          "AEP",
          function(.Object, d, alpha, nb_iteration) {
			if(nargs() < 3) stop("missing arguments")     
			if (is.data.frame(d) == FALSE) stop("the first parameter must be of type 'data.frame'") 			
			.Object@data <- d
			.Object@nb_obs <- dim(d)[1]
            .Object@nb_var <- dim(d)[2]  
			.Object@alpha <- alpha
			if (alpha < 1 / .Object@nb_var || alpha >= 1) stop("'alpha' must be in [1/d, 1)") 			
			.Object@nb_iteration <- nb_iteration
			.Object@pobs_data <- pobs(.Object@data)
			.Object@N <- 2^.Object@nb_var
			.Object@binary_matrix <- array(0, c(.Object@N, .Object@nb_var))			
			.Object@sum_ones <- array(0, c(.Object@N))
			.Object@Hmultiplier <- array(0, c(.Object@N))
			
			EmpiricalCDF <- NULL
			#Compute the ECDF of each variable
			for (i in 1:.Object@nb_var) {  
				.Object@EmpiricalCDF <- c(.Object@EmpiricalCDF, ecdf(.Object@data[,i]))
			}
			.Object@ECDF <- function(i, x) .Object@EmpiricalCDF[[i]](x)

			#Initialize the matrix of binary and the sum of ones in each line
			for (i in 2:.Object@N) {
				.Object@binary_matrix[i,] <- .Object@binary_matrix[i-1,]
				j <- 1
				.Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
				while (.Object@binary_matrix[i, j] > 1) {
					.Object@binary_matrix[i, j] <- 0
					j <- j + 1
					.Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
				}
				.Object@sum_ones[i] <- sum(.Object@binary_matrix[i,])
			}

			#Initialize the multiplier of H()
			for (i in 1:.Object@N) {
				if ((.Object@nb_var  -  .Object@sum_ones[i]) %% 2 == 0) .Object@Hmultiplier[i] <- 1
				else  .Object@Hmultiplier[i] <- -1
			}

			# Initialize the vector m
			.Object@m <- array(0, c(.Object@N))
			.Object@m[1] <- 0
			for (j in 2:.Object@N) {
				if (.Object@sum_ones[j] < 1 / .Object@alpha) {
					if((1 + .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
					else .Object@m[j] <- -1
				}
				else if (.Object@sum_ones[j] == 1 / .Object@alpha) .Object@m[j] <- 0
				else {
					if((.Object@nb_var + 1 - .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
					else .Object@m[j] <- -1
				}
			}
			
            .Object
          }
)
setMethod(f = "show", signature = "AEP",
	definition = function(object) {
		cat("Nb observations =", object@nb_obs, "\nNb variables =", object@nb_var, "\n")
		print(object@data)
	}
)

# Define an HyperCube/Simplex class
setClass("HCSimplex", representation(
	b = "array",
	h = "numeric",
	s = "numeric",
	iteration = "numeric"
	)
)
setMethod("initialize",
          "HCSimplex",
          function(.Object, or, lim, multiplier, it) {
			.Object@b <- or
			.Object@h <- lim
			.Object@s <- multiplier
			.Object@iteration <- it
			.Object
		  }
)



VHmeasure <-function(aep, hc) {
	if (class(aep) != "AEP") 
		stop("The parameter 'aep' is not of class 'AEP'")
	if (length(hc@b) != aep@nb_var) {
		stop("hc@b is too long")		
	}
	samplePoints <- array(0, c(aep@N, aep@nb_var))
	Uval <- array(0, c(aep@N, aep@nb_var))
	for (i in 1:aep@N) {
		samplePoints[i,] <- hc@b + (hc@h * aep@binary_matrix[i,])
		for (j in 1:aep@nb_var) Uval[i,j] <- aep@ECDF(j, samplePoints[i,j])	
	}
	H <- C.n(u=Uval, U=aep@pobs_data)
	Hm <- aep@Hmultiplier * H	
	sum(Hm)
}


Prob <- function(aep, value) {
	nb_of_resolution <- 0
	result <- 0
	spx <- new("HCSimplex", array(0, c(aep@nb_var)), value, 1, 1)
	ListOfSimplexes <- c(spx)	
	while (length(ListOfSimplexes) > 0) {
		spx <- ListOfSimplexes[[1]]
		if (length(ListOfSimplexes) > 1) ListOfSimplexes <- ListOfSimplexes[2:length(ListOfSimplexes)]
		else ListOfSimplexes <- NULL
		hc <- new("HCSimplex", spx@b, aep@alpha * spx@h, spx@s, spx@iteration)
		if (spx@iteration < aep@nb_iteration) {
			new.it <- spx@iteration + 1
			for (j in 2:aep@N) {
				new.b <- spx@b + (aep@alpha * spx@h * aep@binary_matrix[j,])
				new.h <- (1 - (aep@alpha * aep@sum_ones[j])) * spx@h
				new.s <- spx@s * aep@m[j]				
				ListOfSimplexes <- c(ListOfSimplexes, new("HCSimplex", new.b, new.h, new.s, new.it))			
			}
		}
		result <- result + hc@s * VHmeasure(aep, hc)
		nb_of_resolution <- nb_of_resolution + 1
	}
#	print(nb_of_resolution)
	result	
} 

       
			
test.debug <- TRUE
if (test.debug) {
	#test case
	nb.obs <- 30000
	nb.var <- 2
	mydata = matrix(runif(nb.obs * nb.var), c(nb.obs, nb.var))
	mydata <- mydata * 50
	my.aep.test <- new("AEP", data.frame(mydata), alpha=0.75, nb_iteration=8)
	print(Prob(my.aep.test, 25))
	#end of test case
} else {	
	par(mfrow=c(2,2))
	mydata <- read.table("output_10_100.txt", header=TRUE, sep=" ")
	observed.exectime <- array(0, c(dim(mydata)[1]))
	for (i in 1:length(observed.exectime)) observed.exectime[i] <- sum(mydata[i,])
	plot( x = 1:length(observed.exectime), y = observed.exectime, 
		main = "Observed execution time",
		xlab = "#run",
		ylab = "Observed Exec Time",
		type = "h", 
		col = "red")

	observed.ecdf <- ecdf(observed.exectime)
	plot(observed.ecdf,
		main = "Observed ECDF",
		xlab = "#run",
		ylab = "P[ExecTime < X]"
	)

	my.aep <- new("AEP", data.frame(mydata), alpha=0.75, nb_iteration=2)
	my.step <- (max(observed.exectime) - min(observed.exectime)) / 100
	estimated.ecdf <- array(0, c(100))

	index <- 1
	steps <- seq(min(observed.exectime), max(observed.exectime), by=my.step)
	for (i in steps) {
		estimated.ecdf[index] <- Prob(my.aep, i)
		index <- index + 1
	}
	
	plot( x = steps, y = estimated.ecdf, 
		main = "Estimated execution time",
		xlab = "Estimated execution time",
		ylab = "P[ExecTime < X]",
		type = "h", 
		col = "red")
}
